home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / sum / levin_u.c next >
Encoding:
C/C++ Source or Header  |  2001-03-15  |  6.4 KB  |  256 lines
  1. /* sum/levin_u.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <gsl/gsl_math.h>
  22. #include <gsl/gsl_test.h>
  23. #include <gsl/gsl_errno.h>
  24. #include <gsl/gsl_sum.h>
  25.  
  26. int
  27. gsl_sum_levin_u_accel (const double *array, const size_t array_size,
  28.                        gsl_sum_levin_u_workspace * w, 
  29.                        double *sum_accel, double *abserr)
  30. {
  31.   return gsl_sum_levin_u_minmax (array, array_size,
  32.                  0, array_size - 1, w, sum_accel, abserr);
  33. }
  34.  
  35. int
  36. gsl_sum_levin_u_minmax (const double *array, const size_t array_size,
  37.             const size_t min_terms, const size_t max_terms,
  38.             gsl_sum_levin_u_workspace * w,
  39.             double *sum_accel, double *abserr)
  40. {
  41.   if (array_size == 0)
  42.     {
  43.       *sum_accel = 0.0;
  44.       *abserr = 0.0;
  45.       w->sum_plain = 0.0;
  46.       w->terms_used = 0;
  47.       return GSL_SUCCESS;
  48.     }
  49.   else if (array_size == 1)
  50.     {
  51.       *sum_accel = array[0];
  52.       *abserr = 0.0;
  53.       w->sum_plain = array[0];
  54.       w->terms_used = 1;
  55.       return GSL_SUCCESS;
  56.     }
  57.   else
  58.     {
  59.       const double SMALL = 0.01;
  60.       const size_t nmax = GSL_MAX (max_terms, array_size) - 1;
  61.       double noise_n = 0.0, noise_nm1 = 0.0;
  62.       double trunc_n = 0.0, trunc_nm1 = 0.0;
  63.       double actual_trunc_n = 0.0, actual_trunc_nm1 = 0.0;
  64.       double result_n = 0.0, result_nm1 = 0.0;
  65.       double variance = 0;
  66.       size_t n;
  67.       unsigned int i;
  68.       int better = 0;
  69.       int before = 0;
  70.       int converging = 0;
  71.       double least_trunc = GSL_DBL_MAX;
  72.       double least_trunc_noise = GSL_DBL_MAX;
  73.       double least_trunc_result;
  74.  
  75.       /* Calculate specified minimum number of terms.  No convergence
  76.          tests are made, and no truncation information is stored.  */
  77.  
  78.       for (n = 0; n < min_terms; n++)
  79.     {
  80.       const double t = array[n];
  81.       result_nm1 = result_n;
  82.       gsl_sum_levin_u_step (t, n, nmax, w, &result_n);
  83.     }
  84.  
  85.       least_trunc_result = result_n;
  86.  
  87.       variance = 0;
  88.       for (i = 0; i < n; i++)
  89.     {
  90.       double dn = w->dsum[i] * GSL_MACH_EPS * array[i];
  91.       variance += dn * dn;
  92.     }
  93.       noise_n = sqrt (variance);
  94.  
  95.       /* Calculate up to maximum number of terms.  Check truncation
  96.          condition.  */
  97.  
  98.       for (; n <= nmax; n++)
  99.     {
  100.       const double t = array[n];
  101.  
  102.       result_nm1 = result_n;
  103.       gsl_sum_levin_u_step (t, n, nmax, w, &result_n);
  104.  
  105.       /* Compute the truncation error directly */
  106.  
  107.       actual_trunc_nm1 = actual_trunc_n;
  108.       actual_trunc_n = fabs (result_n - result_nm1);
  109.  
  110.       /* Average results to make a more reliable estimate of the
  111.          real truncation error */
  112.  
  113.       trunc_nm1 = trunc_n;
  114.       trunc_n = 0.5 * (actual_trunc_n + actual_trunc_nm1);
  115.  
  116.       noise_nm1 = noise_n;
  117.       variance = 0;
  118.  
  119.       for (i = 0; i <= n; i++)
  120.         {
  121.           double dn = w->dsum[i] * GSL_MACH_EPS * array[i];
  122.           variance += dn * dn;
  123.         }
  124.  
  125.       noise_n = sqrt (variance);
  126.  
  127.       /* Determine if we are in the convergence region.  */
  128.  
  129.       better = (trunc_n < trunc_nm1 || trunc_n < SMALL * fabs (result_n));
  130.       converging = converging || (better && before);
  131.       before = better;
  132.  
  133.       if (converging)
  134.         {
  135.           if (trunc_n < least_trunc)
  136.         {
  137.           /* Found a low truncation point in the convergence
  138.              region. Save it. */
  139.  
  140.           least_trunc_result = result_n;
  141.           least_trunc = trunc_n;
  142.           least_trunc_noise = noise_n;
  143.         }
  144.  
  145.           if (noise_n > trunc_n / 3.0)
  146.         break;
  147.  
  148.           if (trunc_n < 10.0 * GSL_MACH_EPS * fabs (result_n))
  149.         break;
  150.         }
  151.  
  152.     }
  153.  
  154.       if (converging)
  155.     {
  156.       /* Stopped in the convergence region.  Return result and
  157.          error estimate.  */
  158.  
  159.       *sum_accel = least_trunc_result;
  160.       *abserr = GSL_MAX_DBL (least_trunc, least_trunc_noise);
  161.       w->terms_used = n;
  162.       return GSL_SUCCESS;
  163.     }
  164.       else
  165.     {
  166.       /* Never reached the convergence region.  Use the last
  167.          calculated values.  */
  168.  
  169.       *sum_accel = result_n;
  170.       *abserr = GSL_MAX_DBL (trunc_n, noise_n);
  171.       w->terms_used = n;
  172.       return GSL_SUCCESS;
  173.     }
  174.     }
  175. }
  176.  
  177.  
  178. int
  179. gsl_sum_levin_u_step (const double term, const size_t n, const size_t nmax,
  180.               gsl_sum_levin_u_workspace * w, double *sum_accel)
  181. {
  182.  
  183. #define I(i,j) ((i)*(nmax+1) + (j))
  184.  
  185.   if (n == 0)
  186.     {
  187.       *sum_accel = term;
  188.       w->sum_plain = term;
  189.  
  190.       w->q_den[0] = 1.0 / term;
  191.       w->q_num[0] = 1.0;
  192.  
  193.       w->dq_den[I (0, 0)] = -1.0 / (term * term);
  194.       w->dq_num[I (0, 0)] = 0.0;
  195.  
  196.       w->dsum[0] = 1.0;
  197.  
  198.       return GSL_SUCCESS;
  199.     }
  200.   else
  201.     {
  202.       double result;
  203.       double factor = 1.0;
  204.       double ratio = (double) n / (n + 1.0);
  205.       unsigned int i;
  206.       int j;
  207.  
  208.       w->sum_plain += term;
  209.  
  210.       w->q_den[n] = 1.0 / (term * (n + 1.0) * (n + 1.0));
  211.       w->q_num[n] = w->sum_plain * w->q_den[n];
  212.  
  213.       for (i = 0; i < n; i++)
  214.     {
  215.       w->dq_den[I (i, n)] = 0;
  216.       w->dq_num[I (i, n)] = w->q_den[n];
  217.     }
  218.  
  219.       w->dq_den[I (n, n)] = -w->q_den[n] / term;
  220.       w->dq_num[I (n, n)] =
  221.     w->q_den[n] + w->sum_plain * (w->dq_den[I (n, n)]);
  222.  
  223.       for (j = n - 1; j >= 0; j--)
  224.     {
  225.       double c = factor * (j + 1) / (n + 1);
  226.       factor *= ratio;
  227.       w->q_den[j] = w->q_den[j + 1] - c * w->q_den[j];
  228.       w->q_num[j] = w->q_num[j + 1] - c * w->q_num[j];
  229.  
  230.       for (i = 0; i < n; i++)
  231.         {
  232.           w->dq_den[I (i, j)] =
  233.         w->dq_den[I (i, j + 1)] - c * w->dq_den[I (i, j)];
  234.           w->dq_num[I (i, j)] =
  235.         w->dq_num[I (i, j + 1)] - c * w->dq_num[I (i, j)];
  236.         }
  237.  
  238.       w->dq_den[I (n, j)] = w->dq_den[I (n, j + 1)];
  239.       w->dq_num[I (n, j)] = w->dq_num[I (n, j + 1)];
  240.     }
  241.  
  242.       result = w->q_num[0] / w->q_den[0];
  243.  
  244.       *sum_accel = result;
  245.  
  246.       for (i = 0; i <= n; i++)
  247.     {
  248.       w->dsum[i] =
  249.         (w->dq_num[I (i, 0)] -
  250.          result * w->dq_den[I (i, 0)]) / w->q_den[0];
  251.     }
  252.  
  253.       return GSL_SUCCESS;
  254.     }
  255. }
  256.